Setup
Data describing commercial landings (tonnes) in 2017.
# Read in commercial landings data
fish <- rast(here("Reference", "CRS", "commercial_landings_2017.tif"))
## Take a look
plot(fish)How many fish were captured?
## sum
## commercial_landings_2017 86331119
Check against Sea Around Us to make sure this seems generally reasonable.
What is the area of each cell?
Describe what the plotted image is saying?
Given this, what is kind of deceptive about these plots?
cellSize()Compute the area covered by individual raster cells. Computing the surface area of raster cells is especially relevant for longitude/latitude rasters.
But note that for both angular (longitude/latitude) and for planar (projected) coordinate reference systems raster cells sizes are generally not constant, unless you are using an equal-area coordinate reference system.
For planar CRSs, the area is therefore not computed based on the linear units of the coordinate reference system, but on the actual area by transforming cells to longitude/latitude. If you do not want that correction, you can use transform=FALSE or init(x, prod(res(x)))
Let’s convert these data to another coordinate reference system:
moll <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
fish_moll <- terra::project(fish, crs(moll))
# take a look
par(mar = c(5, 4, 4, 0)) # set margins: bottom, left, top, right
plot(fish_moll)## sum
## commercial_landings_2017 389271065
# 389.3 million
# calculate cell area
cell_area_moll <- cellSize(fish_moll)
#par(mai = c(1, 1, 3, 5)) # set margins: bottom, left, top, right
#par(oma = c(1, 1, 3, 4))
#par(mar = c(1,1,1,1))
plot(cell_area_moll, main = "Mollweide projection cell size")## class : SpatRaster
## dimensions : 764, 1546, 1 (nrow, ncol, nlyr)
## resolution : 23330.19, 23330.19 (x, y)
## extent : -18037447, 18031019, -8861637, 8962624 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source(s) : memory
## name : area
## min value : 0
## max value : 547291328
## class : SpatRaster
## dimensions : 347, 720, 1 (nrow, ncol, nlyr)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -85.5, 88 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : commercial_landings_2017
## name : area
## min value : 122442386
## max value : 3077249667
Compare to the original tonnes? What is happening? How can you improve this? - original tons: 86331119 - Mollweide tons: 389271065 - the cell size area is now relatively consistent across the globe – old cell sizes have stretched etc., to fit a standardized cell size (I tried to extend the margins to see the full scalebar label values but struggled) - perform global calculations before reprojecting
could take raster of cell size, multiply by original (if density) raster, then sum values
^ divide if count
tonnes / area = density, then project
take new cell size, multiply by density raster
sum at end to check again – might need to reduce resolution, summary before projecting data
# reproject to Mollweide
fish_density_moll <- terra::project(fish_density, crs(moll))
fish_density_moll## class : SpatRaster
## dimensions : 764, 1546, 1 (nrow, ncol, nlyr)
## resolution : 23330.19, 23330.19 (x, y)
## extent : -18037447, 18031019, -8861637, 8962624 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source(s) : memory
## name : commercial_landings_2017
## min value : 0.000000e+00
## max value : 4.908811e-05
# multiply density by area to get count (tonnes)
new_tonnes <- fish_density_moll * new_cell_size
plot(log(new_tonnes + 1))## sum
## commercial_landings_2017 85253920
## commercial_landings_2017
## Min. : 0.00
## 1st Qu.: 2.43
## Median : 22.25
## Mean : 651.53
## 3rd Qu.: 84.65
## Max. :125112.30
## NA's :47150
# Get the cell numbers
cells <- 1:terra::ncell(fish)
# Get the x and y coordinates for each cell
xy <- terra::xyFromCell(fish, cells)
# Extract the values
values <- terra::values(fish)
# Combine everything into a data frame
fish_df <- data.frame(
x = xy[,1],
y = xy[,2],
value = values
)
fish_no_na <- fish_df %>%
na.omit()See what happens when we make multiple CRS conversions:
# Robinson projection
rob <- "+proj=robin +datum=WGS84 +units=m +no_defs"
fish_moll_rob <- terra::project(fish_moll, crs(rob))
plot(log(fish_moll_rob + 1))# now back to lat long:
fish_moll_rob_latlon <- project(fish_moll_rob, fish)
plot(log(fish_moll_rob_latlon + 1))# check to see if it exactly matches the start
check <- fish_moll_rob_latlon - fish
plot(log(check + 1))## class : SpatRaster
## dimensions : 347, 720, 1 (nrow, ncol, nlyr)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -85.5, 88 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : commercial_landings_2017
## name : commercial_landings_2017
## min value : -66229.40
## max value : 35487.52
What is this telling us?
the values have changed notably between projections
some are highly inflated (large positive numbers), other deflated (large negative in new projection)
adding error when we flip between resolutions
they’re all just approximations (by nature)
## class : SpatRaster
## dimensions : 772, 1542, 1 (nrow, ncol, nlyr)
## resolution : 22047.33, 22047.33 (x, y)
## extent : -17003404, 16993576, -8462847, 8557690 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source(s) : memory
## name : commercial_landings_2017
## min value : 0.0
## max value : 132159.3
fish_moll_rob_latlon <- project(fish_moll_rob, fish) # match template (convert crs and resolution)
fish_moll_rob_latlon## class : SpatRaster
## dimensions : 347, 720, 1 (nrow, ncol, nlyr)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -85.5, 88 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : commercial_landings_2017
## name : commercial_landings_2017
## min value : 0.0
## max value : 104592.5
fish_moll_rob_latlon_v2 <- project(fish_moll_rob, crs(fish)) # match crs (just convert crs)
fish_moll_rob_latlon_v2## class : SpatRaster
## dimensions : 750, 1553, 1 (nrow, ncol, nlyr)
## resolution : 0.2316815, 0.2318308 (x, y)
## extent : -180, 179.8014, -85.87307, 88 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : commercial_landings_2017
## min value : 0.0
## max value : 127646.3